{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RANSAC算法的核心思想\n",
    "它的核心思想是不断随机采样少量点去拟合一个模型。然后用其他数据点判断这个模型是否足够好。判断这个模型是否足够好的标准就是有多少个点满足这个模型。如果这个模型足够好的话那么用满足这个模型的那些点去重新拟合一个新的模型。（不满足条件的那些点就扔了。这个和最小二乘法不同，最小二乘法是把所有点都放一起取平均）\n",
    "\n",
    "基本概念：每次最少采样点数m（一条直线最少需要m=2个点），判断其他数据点是否满足模型时候用的最大容忍误差e（当代入到m个点所拟合的那个模型中发现误差大于e就认为是不合格的点outliers），足够好的模型的最少合格点数n（如果满足这个模型的合格点大于等于n就认为这个模型足够好了终止迭代），迭代次数K(如果一直找不到足够好的模型则不断迭代采样直到迭代次数截止为止)。\n",
    "\n",
    "一般：$n-m>5，K=2p^{(-m)}$,其中p是指随机的一个数据点落在拟合模型范围内的概率（这个概率是由最大容忍误差决定的），当然K你随便取也可以。\n",
    "下面是Python代码实践用RANSAC去除离群点然后拟合一条曲线："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAGs1JREFUeJzt3W+MXGd1x/Hv8XrA40CyQVlovM5ioyKnbUJtOkpRV6KQAE5JlFipSqANQoXKQgIa3GKwiyrgTbOqWxKkVqqskAIiAiNiTESgTloHUSISWGODCYkRIiF4HWqjZPkjD7Ben77YHWd29t479+69M3fmub+PZMU7e/fOMxP5zLPnOc95zN0REZFwrCp7ACIiUiwFdhGRwCiwi4gERoFdRCQwCuwiIoFRYBcRCYwCu4hIYBTYRUQCo8AuIhKY1WU86SWXXOIbNmwo46lFRIbW4cOHf+buY92uKyWwb9iwgenp6TKeWkRkaJnZj9Ncp1SMiEhgFNhFRAKjwC4iEhgFdhGRwCiwi4gEppSqGBGRXjlwZIY9B49zcrbJutE6O7duYtuW8bKH1VcK7CISjANHZti9/xjNuXkAZmab7N5/DKBSwT11KsbM7jKzU2b2vYjvvc/M3MwuKXZ4IiLp7Tl4/HxQb2nOzbPn4PGSRlSOLDn2TwDXdj5oZpcBrweeKmhMIiIrcnK2menxUKUO7O7+NeCZiG/dDrwf0KnYIlKqdaP1TI+HKldVjJndAMy4+3cKGo+IyIrt3LqJem1kyWP12gg7t24qaUTlWPHiqZmtBT4IvCHl9duB7QATExMrfVoRkVitBdKqV8WYe/oMipltAL7k7leY2ZXA/wBnFr+9HjgJXOXuP026T6PRcDUBExHJxswOu3uj23UrnrG7+zHgxW1P+CTQcPefrfSeIiKSX5Zyx88A3wA2mdkJM3tH74YlIiIrlXrG7u5v6fL9DblHIyIiualXjIhIYBTYRUQCo8AuIhIYBXYRkcAosIuIBEaBXUQkMArsIiKBUWAXEQmMAruISGAU2EVEAqPALiISGAV2EZHAKLCLiARGgV1EJDAK7CIigVFgFxEJTJYTlO4ys1Nm9r22x/aY2eNm9l0z+4KZjfZmmCIiklaWGfsngGs7HnsAuMLdXwH8ANhd0LhEpGIOHJlhcuoQG3fdx+TUIQ4cmSl7SEMrdWB3968Bz3Q8dr+7n1388mFgfYFjE5GKOHBkht37jzEz28SBmdkmu/cfU3BfoSJz7G8HvhL3TTPbbmbTZjZ9+vTpAp9WRIbdnoPHac7NL3msOTfPnoPHSxrRcCsksJvZB4GzwN1x17j7XndvuHtjbGysiKcVkUCcnG1melyS5Q7sZvY24Hrgr9zd8w9JRKpm3Wg90+OSLFdgN7NrgQ8AN7j7mWKGJCJVs3PrJuq1kSWP1Wsj7Ny6qaQRDbfVaS80s88ArwEuMbMTwIdYqIJ5PvCAmQE87O7v7ME4RSRg27aMAwu59pOzTdaN1tm5ddP5x4fdgSMzfX1tVkb2pNFo+PT0dN+fV0Sk31oVP+2Lw/XaCLfddGXm4G5mh9290e067TwVEemhMip+FNhFRHqojIofBXYRkR4qo+JHgV1EpIfKqPhJXRUjIhKKflaplFHxo8AuIpXSWaXS6ksD9DS497N0U6kYEamUKvSlUWAXkUqpQl8aBXYRqZQq9KVRYBeRSqlCXxotnopIpfSrSqXf/WHaKbCLSOX0ukqljMqbdgrsIiIxVjrrTqq8UWAXESlJnll32ZU3CuwiIm1as/SZiCCcdta9brQe+fP9qrxRVYyIyKLWLD0qKLekmXWXXXmT5QSlu1g42/SUu1+x+NiLgH3ABuBJ4E3u/mzxwxQR6b2o3HgnBzZ/5H7MYPbMXGTuvewToVKfoGRmrwZ+BXyqLbD/M/CMu0+Z2S7gYnf/QLd76QQlERlEG3fdx0rOlFvpiUhZFX6Ckrt/DXim4+EbgU8u/v2TwLbUIxQRGQAHjswwOXWIjbvuY9XC2c2ZDVqvmbyLpy9x96cB3P1pM3tx3IVmth3YDjAxMZHzaUVE8uusfJmPyGDUayNd0zOwkHsvc1NSu74tnrr7XndvuHtjbGysX08rIhIrLqc+YoYB46N1brvpSsZTVLM4sGPfUWZmmzjPlUceODJT+Li7yTtj/z8zu3Rxtn4pcKqIQYmI9ENchcs5d56Yum7JY+0z+zid8/1+bkpql3fGfi/wtsW/vw34Ys77iYj0TdpOj9u2jJ+fuRswWq9x8dpaqucoox1wlnLHzwCvAS4xsxPAh4Ap4HNm9g7gKeAvejFIEZFe2Ll107KZeFy9eVR/mTRVNGW0A04d2N39LTHfuqagsYiI9FXeevO4HaYtZbUDVksBEam0PJ0eo2b8xkKufbzEqhgFdhEZemnKDHtRilj2DtM4CuwiMtTSdGHsZX/0Xvd2Xwk1ARORoZbU+zzLNSHRjF1EhlpS7/OkFrxJPzvsNGMXkaEWV054Ub3WtQVvGaWI/aDALiIDrb1J1+TUoWVb9ON6n5uRuFO0rFLEflBgF5GB1X7wRVz/lc5doa3+LrNn5mLv27pm0BY9i6Icu4j0XC8OhW59P+6ecbn18dE6D+26OucrGmyasYtIT6WZdceJW9xs3SPpnmUfT1cmBXYR6ak8pYZxi5sjZl3vGZeiCTX90k6pGBHpqaRyxDjtZYqtLfotSQdfdN5zEDcP9YMCu4gUJiqXHtcoK2423rlL1FnefyUufx5q+WJWSsWISCHicumvvXwsU647KnXTCuoP7bqabVvGK50/T0Mz9oANyvmLw0TvWTbt79cqs2Vnhjbn5nnw8dPcdtOV56+7qF7DbOEYuQ/f+yhmMHtm7vz7nSZ1M6jNtwaFecThrZlvYrYD+BsWPliPAX/t7r+Ou77RaPj09HTu55V4nb/OwsKMpiqLRyuh9yybqPcrisH5Y+a6/Uy9NsKa2iqejahBr0KZYjdmdtjdG92uy52KMbNx4G+BhrtfAYwAb857X8mnak2PiqD3LJu4g6A7tee9u/1Mc24ed5RmyamoHPtqoG5mq4G1wMmC7isrtJJKhKrTe7ZUt638ad4XYyHX3vr5ND/z8+ZcZcsUi5I7x+7uM2b2LyycedoE7nf3+3OPTHLJWokges/apelfHvd+jSzm2tvLFFs/P7q2FplmabdutF7ZMsWiFJGKuRi4EdgIrAMuMLNbIq7bbmbTZjZ9+vTpvE8rXahqIDu9Z8+JS0u9d9/R87PvuPfrX9/0h4yP1pcd8hyXZun8+Sq+30UrIhXzOuAJdz/t7nPAfuBPOi9y973u3nD3xtjYWAFPK0mqvOtupfSePScpZTIz22THvqO8d99Rnr96FRevrS17v+J+vjPNMlqvRf685FNEueNTwKvMbC0LqZhrAJW8DAD9OptdaO/ZSss349IsLa3Z+GxzjnpthNtv3rzkvklprdDe40GUe8bu7o8Anwe+zUKp4ypgb977ikg+eZpvRaVZ4kRVDimtVa5CNii5+4eADxVxLxEpRlL5ZvuMOWlWn3SsXLuoHi2tn9cGov4rZINSVtqgJNI73c75hOd6rgDLNgx19mWJuibqflXfPNQPaTcoKbAHRlviqy3tblBI3uXZfs1tN10JkNhtUYue/ZE2sKtXTEDS1B5L2NLuBoWFtEy3a1upm1bzLSh38qCJSzoK7AFJm1OVMEQFuV7skh2UHueauKSnwB4QbYkPU1QAByKDXNzOzpGIzouwUEf+m7PnEmfug7LzVhOX9BTYA6It8eGIO0GoFcDX1FZFBrnm3HxkDvzP/2icew7PLOtc+eEb/gBIzp8PSomiJi7pKbAHZOfWTZFtZwflH6akE3WCULtuufGoE4e2bRmn8dIXxeanByF/3o0mLumpKiYwg/wPU5aK+381OXUoVe14N6GVIKpfvqpiKqsfC1v68MgvaSEwTWohTW48tBSFNj2lp8AumagyIbuoD8KkhcBufVqicuNRQkxRqM9MOjrMWjLRKUPZxPVriQvGJ2ebkX1WbPG/o/Uaa2qr2LHvKHsOHmfn1k3ccfNm9WWRJRTYJRNVJmQT90E4YhZ5fav7YWf74Ntv3swdN2/mN2fP8eyZuSUfEoDaDcsSSsVIJqpMyCbuA2/enXptJLaCKSrlMDl1KPa3pfadoSKasUsmaseaTdIHXtwhFXH025KkpcAumeiUoWyS+prPNuf49dw5br95c6oZd9yHhH5bkk5KxUhmVa9MyFLu2a2vefuW+G731QY0SauQwG5mo8CdwBUsbHh7u7t/o4h7iwySpHJPiK6xbv3ZuOu+ZbtIYSGVkqaMVHXcklYhO0/N7JPA/7r7nWb2PGCtu8/GXa+dpzIM4urPo2beURuGOndFxu0oHV9MpcR9L6Tdo5JP3w7aMLMLge8AL/OUN1Ngl0EXt309ba/zdkknEbWC/459RyNn8wY8MXVd5ueUMKUN7EUsnr4MOA38p5kdMbM7zeyCAu4r0jMHjswwOXWIjbvuY3Lq0LIDnrPWnydJU2+uhVEpUhE59tXAK4H3uPsjZvYxYBfwj+0Xmdl2YDvAxMREAU8rkk23VrjwXB47S/15Gt3qzbUwKkUqYsZ+Ajjh7o8sfv15FgL9Eu6+190b7t4YGxsr4GlF0mvf2g/RrXDb2yLEzZRbs+zxFcykk+rNVUYqRco9Y3f3n5rZT8xsk7sfB64Bvp9/aCLFSXMWaHvgTZpBt6pcsrbX7ZZWqXoZqRSnqA1K7wHuNrPvApuBfyroviKFSLM7sz3wpplBx+3CveVVE9qdK6UqpI7d3Y8CXVdqyxRyD/GQX1tR0rTC7Qy87TPo1nu8Y9/RZe9x1HufdFqRSK9V4gSlkE9eCfm1FSnqfYo6Pi7tz+o9ljLoBKU2IZ9uHvJrK1KeXZt6j2XYVCKwh9wVL+TXVrSVLk7qPZZhU4nAHnIP8ZBfW5LOdYXXXj7Gg4+f7klOu6rvsQyvSrTtDbmHeMivLU7UcXOffvipZcfPde4mTbpf0i7UKr7HMtwqMWMPuSteyK8tTpqa9LQ5cHVVlBBVoipGBl+Wks249red0jTQSuq4qK6KMmj62QRMJJeo1EpSKiVtbjvNdVoYlRApsEvpksoJoyQdN9eSNgeurooSIgV26auohcqss+ao7f63vGri/Nej9RpraqvYse9o5GJoOy2MSogqsXgq5erWLnd0bY1nz8wt+7mkWXNcTXqaxdDO+4AWRiUsWjyVnorajt8p6Vg5yBZ0tRgqIdPiqQyENKWJP2/ORXZSBJYtqu7Yd5QNMfXmoMVQEVAqRnLqVqaYtl1uVGplcurQsg+FpFOPWvfSLlGpOs3YZcXSlCl2C6hRC5WtBdZuh1hEVc5oMVREM3ZJKWpmnqbrYdRJREntctPk5Nt1/kagxVCRAgO7mY0A08CMu19f1H2lfHGVJnHBtz3YZg20aXLy7aJ+I9ARc1J1Rc7YbwUeAy4s8J4yAOJm5iNmzEdUVa0bra/4VKeknHx7qSQoxSISp5Acu5mtB64D7izifjJY4oLtvHtkPvu1l49lahHQLi4nPz5a5/abNyeeQSoiC4paPL0DeD9wLu4CM9tuZtNmNn369OmCnlb6ISnYRpUpPvj46UwtAtolLX5u2zLOQ7uu5omp63ho19UK6iIxcqdizOx64JS7Hzaz18Rd5+57gb2wsEEp7/NK/0QtgLYH284Au2Pf0cj7pCl91OKnSH5F5NgngRvM7I3AGuBCM/u0u99SwL1lAGQNtnlrybX4KZJP7sDu7ruB3QCLM/b3KaiHIesCaFxPGNBCp0g/qY5dImVtptV5vZNcry7lWmnVkgyHQgO7u38V+GqR95T+ap91d0o6bi6qJLIV1NV8a7Bk/dCW4aOWAnJee4uAOFmbbKn51uDJerCJDB+lYiqu/VfyVTEbjtolnTik5lvDQR/C4dOMvcI6m3h1C+pJC6BqvjU8dBxg+BTYKyxLX5ZuOz2jjqvTztDBpA/h8CkVU2FpfvVunWSUJkCr/nw4aBNY+BTYK6iVV49LvIyYcc5d/+ADpg/hsCmwD5Eiao+79TvPMkMXkcGkw6wHRLegHRWQV7IBKOlkohA3EmkjjoQk7WHWmrEPgDQbRuI2AMVdHycur24Q3EYibcSRqlJVzABIs2Gk20Jn2g0mVSp100YcqSoF9gEQF7RnZptMTh3iwJGZVIE3TZVLlUrdtBFHqkqBvQ8OHJlhcuoQG3fddz5Qt0sK2q30wWsvH1sWkDulCf5Vqjev0m8nIu2UY++xNHneqIMs2jXn5nnw8dPcdtOVhbTFrUqpW9IBISIhU2AvUFQFRlKetxVc2xdI4ypWTs42lwRkVXt0p404UlUqdyxIVDlivTaSuGU/qrwwrhxR7W9FJG25Y+4cu5ldZmYPmtljZvaomd2a957DKG5mPmIW+zMzs0127DvKhrbce5UWN0WkN4pYPD0L/L27/x7wKuBdZvb7Bdx3qMRVWsy7Jy56RtWiV2VxU0R6o4gzT58Gnl78+y/N7DFgHPh+3nsPurS9zJ+/ehVraqt49sxc4v1aufeHdl2tQC4iK1bo4qmZbQC2AI8Ued9B1JlTT+plPtuco14b4eK1ta7BXTXWIpJXYYHdzF4A3AO8191/EfH97cB2gImJiaKetu+SzgSFhc6IUUG+OTfP81ev6rqgqhprEcmrkA1KZlZjIajf7e77o65x973u3nD3xtjYWBFPW4hum4c6r+12Jug5d+KWS3/enDufPweWXadFUhEpQu4Zu5kZ8HHgMXf/aP4h9U/WJlFpThxqzbjjzv9ULbqI9FoRqZhJ4K3AMTM7uvjYP7j7lwu4d0+l2TwE3dMvLe0z7jQ7HquyA1RE+quIqpivszyrMNC6Ber2BcxuB1O0RG020mxcRMpQuZYCaQJ1+wJmt/RL3IlDmo2LSFkqF9jTBOr2lElS+WGIJw6JyPCrXGDPGqjXjdbVu0VEhkrl+rHH1Ym3AnXn7Fu9W0Rk2FRuxp7Uo7u9/PCieg0zmD0zx0X1Gmtqq5g9M6eFUBEZeJUL7HE9umFpieJs87mt/62WALffvFkBvUCq4xfpjSACe9YAEVWxMjl1KHFRNaq+XVYu6+YwEUlvaAN7ey16+zFxKw0QaZpvqUFXcdJuDhOR7IZy8bSzZ0tny61WgMgiTfMtNegqTtyHpD48RfIbysCepmdL1gARVf3STpUwxYr7kNSHp0h+QxnY0wTtrAFi25bxJScXjdZrXLy2FswpRlm6WPaDykhFemcoc+xxm4ZaosoXV7qoGoJBXKiMq04K8f0X6TfzhJN/eqXRaPj09PSKfz6q30trAXU8pnwR4vu6hG5y6pB2z4oEwMwOu3uj23VDOWNPmu0ldW6satWFFipFqmUoAztEp03SdG6sYjCLS11poVIkTEO5eBonywlHVaKFSpFqKWTGbmbXAh8DRoA73X2qiPt26rYY2m02XtVgpoVKkWop4szTEeDfgdcDJ4Bvmdm97v79vPdul6ayI6lapuq900Ot+BGR5YpIxVwF/NDdf+TuvwU+C9xYwH2XSNqC3hKXcrjj5s2RLXlFREJURCpmHPhJ29cngD8u4L5LpKnsUMpBRKSYwB51kPWy4ngz2w5sB5iYmMj8JGkrO5RyEJGqKyIVcwK4rO3r9cDJzovcfa+7N9y9MTY2lvlJVNkhIpJOETP2bwEvN7ONwAzwZuAvC7jvEkqziIikkzuwu/tZM3s3cJCFcse73P3R3COLoDSLiEh3hdSxu/uXgS8XcS8REcknqJ2nIiKiwC4iEhwFdhGRwCiwi4gERoFdRCQwCuwiIoFRYBcRCYwCu4hIYBTYRUQCo8AuIhIYBXYRkcAosIuIBEaBXUQkMArsIiKBKaRtr/TOgSMzOlxERDJRYB9gB47MsHv/MZpz8wDMzDbZvf8YgIK7iMTKlYoxsz1m9riZfdfMvmBmo0UNTBaOAWwF9Zbm3Dx7Dh4vaUQiMgzy5tgfAK5w91cAPwB25x+StJycbWZ6XEQEcgZ2d7/f3c8ufvkwsD7/kKRl3Wg90+MiIlBsVczbga8UeL/K27l1E/XayJLH6rURdm7dVNKIRGQYdF08NbP/Bn4n4lsfdPcvLl7zQeAscHfCfbYD2wEmJiZWNNiqaS2QqipGRLIwd893A7O3Ae8ErnH3M2l+ptFo+PT0dK7nFRGpGjM77O6NbtflKnc0s2uBDwB/mjaoi4hIb+XNsf8b8ELgATM7amb/UcCYREQkh1wzdnf/3aIGIiIixVCvGBGRwCiwi4gEJndVzIqe1Ow08OMct7gE+FlBwxkGVXu9oNdcBVV7vZD/Nb/U3ce6XVRKYM/LzKbTlPyEomqvF/Saq6Bqrxf695qVihERCYwCu4hIYIY1sO8tewB9VrXXC3rNVVC11wt9es1DmWMXEZF4wzpjFxGRGEMV2M3sWjM7bmY/NLNdZY+n18zsLjM7ZWbfK3ss/WJml5nZg2b2mJk9ama3lj2mXjKzNWb2TTP7zuLr/UjZY+oXMxsxsyNm9qWyx9IPZvakmR1bbL/S0y6IQ5OKMbMRFk5pej1wAvgW8BZ3/36pA+shM3s18CvgU+5+Rdnj6QczuxS41N2/bWYvBA4D20L9/2xmBlzg7r8ysxrwdeBWd3+45KH1nJn9HdAALnT368seT6+Z2ZNAw917Xrs/TDP2q4AfuvuP3P23wGeBG0seU0+5+9eAZ8oeRz+5+9Pu/u3Fv/8SeAwItgG9L/jV4pe1xT/DMdvKwczWA9cBd5Y9lhANU2AfB37S9vUJAv4HL2BmG4AtwCPljqS3FlMSR4FTwAPuHvTrXXQH8H7gXNkD6SMH7jezw4sHD/XMMAV2i3gs+JlNVZnZC4B7gPe6+y/KHk8vufu8u29m4czgq8ws6LSbmV0PnHL3w2WPpc8m3f2VwJ8B71pMtfbEMAX2E8BlbV+vB06WNBbpocVc8z3A3e6+v+zx9Iu7zwJfBa4teSi9NgncsJhz/ixwtZl9utwh9Z67n1z87yngCyykl3timAL7t4CXm9lGM3se8Gbg3pLHJAVbXEz8OPCYu3+07PH0mpmNmdno4t/rwOuAx8sdVW+5+253X+/uG1j4d3zI3W8peVg9ZWYXLBYDYGYXAG8AelbtNjSB3d3PAu8GDrKwoPY5d3+03FH1lpl9BvgGsMnMTpjZO8oeUx9MAm9lYRZ3dPHPG8seVA9dCjxoZt9lYfLygLtXovyvYl4CfN3MvgN8E7jP3f+rV082NOWOIiKSztDM2EVEJB0FdhGRwCiwi4gERoFdRCQwCuwiIoFRYBcRCYwCu4hIYBTYRUQC8//HIosXi6FougAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def generate_data():\n",
    "    x = np.linspace(0,5,100)\n",
    "    y = 2*x+np.random.rand(x.shape[0])\n",
    "    # 制造离群点\n",
    "    outlier_num = 20\n",
    "    for i in range(outlier_num):\n",
    "        y[np.random.randint(100)] += 10*(np.random.rand()-0.5)\n",
    "    return x,y\n",
    "x,y  = generate_data()\n",
    "plt.scatter(x,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAFehJREFUeJzt3W+MXOV1x/Hf8TKpByhsKtwUBlxTKTJtcIvDKqJdNQpLWtOAiOWoJZESpVUkv0lbIC3R8or0ReWVqBL6qpKVpE1FBE6BuhFIcaKYKsIqJGvslIChisKfsNB6o3iTEK9gMacvdsYeX987c+f+//P9SJbt8fx5RgnHx+c5z3nM3QUAqL8NZS8AAJANAjoANAQBHQAagoAOAA1BQAeAhiCgA0BDENABoCEI6ADQEAR0AGiI84r8sEsuucS3bNlS5EcCQO0dPnz4J+6+adzzCg3oW7Zs0eLiYpEfCQC1Z2YvxXkeJRcAaAgCOgA0BAEdABqCgA4ADUFAB4CGIKADQEMU2rYIAG2y/8iS7jnwvF5dWdVl013duWOrdm7v5fZ5BHQAyMH+I0u66+Gntbp2SpK0tLKqux5+WpJyC+qUXAAgB/cceP50MB9YXTulew48n9tnEtABIAevrqxO9HgWKLkAQIQ0NfDLprtaCgnel013s17maWToABBiUANfWlmV60wNfP+RpVivv3PHVnU7U2c91u1M6c4dW3NY7ToCOgCESFsD37m9pz27tqk33ZVJ6k13tWfXNrpcAKAogzJLWLlEmqwGvnN7L9cAHkRAB4C+YKthmDxr4GkR0AGgL6zMMmxQAy/6wFBcYwO6mX1Z0s2Sjrv71f3Hfk3SPklbJL0o6c/c/UR+ywSA/I0qp0x3OzKTbt93VCbJ+48XcWAorjibov8i6cbAY/OSvu3u75b07f7vAaDWosop092O3njrbZ04uSbpTDAfyPvAUFxjA7q7f0fSTwMPf1jSV/q//oqknRmvCwAKF9VqaKaRpRgp3wNDcSVtW3yXu78mSf2ffz3qiWa228wWzWxxeXk54ccBQP6iWg1X+pn5KFXYLM19U9Td90raK0kzMzPBf6kAQOnGbXKOamOUzj4wVOaGadIM/f/M7FJJ6v98PLslAUBx4pwIDSvFWP/n4QNDaU+XppU0oH9d0if7v/6kpP/IZjkAUKw4J0LDSjFfuPUavbhwkw7Nz0mSZhcO6vZ9RwufsDgsTtvi/ZI+IOkSM3tF0t2SFiR9zcw+JellSX+a5yIBIC9xpyJGnfqMcxipqA3TsQHd3T8W8Uc3ZLwWAChc2qmI4w4jTfJeaTGcC0CrpZ2KOC77znvC4jCO/gNotUEZJeu559J6rb3ILhcCOoDWSzMV8c4dW8+poXc7U7mPyg1DQAeAFNJm+FkioANASkXPPY/CpigANAQZOoDaq+p88qIR0AHUWvBgz6j55E0P/JRcANRa3Mucy56zUgQCOoBai3t0P27grzMCOoBaizpWH3w8buCvMwI6gFqLe3Q/buCvMwI6gFqLumUouNmZdmZLHdDlAqD2Rh3sGe5subjb0cbOBq2cXGtklwsBHUBjBVsaV1bX1O1M6Qu3XtOoQD5AyQVAY7Whs2UYGTqASktzGKgNnS3DyNABVFbaw0Bt6GwZRkAHkKv9R5Y0u3BQV84/qtmFgxOdzExbMmlDZ8swSi4AcjPJnJUwo0omcUoxVZpVXgRz98I+bGZmxhcXFwv7PADlml04GHo9W2+6q0Pzc4lfP93t6I233q7ELUFFMLPD7j4z7nmUXADkJu2mZFTJxEyt6l6Ji4AOIDdpNyWjToGunFwLfX5Tu1fiooYOIDdRFyiP25QcVx+/58DzoaWYpnavxEWGDiA3ceesDIvTqti27pW4yNABZCIqq570AuVRrYqD92lb90pcBHQAqaVtTxwWdyN10r8o2oCSC4DUspyZ0rbTnVkiQweQWtysOliWuf6qTXrsueWzyiZJN1JBhg4ghcGx/qjjicNZddhm531PvHzO5qekiTdSsY4MHUAiwbp5UDCrDivLBA3KNIfm5wjgCaTK0M3sDjN7xsx+YGb3m9nGrBYGoNpGBeiwrDruoZ+2Hw5KI3FAN7OepL+WNOPuV0uakvTRrBYGoNqiAq9JoRl23E1NNj+TS1tyOU9S18zWJJ0v6dX0SwJQBeNOa1423Y08rRn22rDNziA2P9NJnKG7+5Kkf5D0sqTXJP3M3b+Z1cIAlCfNac3rr9oU+lrp3M3Oj1+3mc3PDCUen2tm75T0kKRbJa1I+jdJD7r7fYHn7Za0W5I2b9587UsvvZRqwQDyF3fsbVgmHjVnJe7IXJwr7vjcNCWXD0p6wd2X+x/4sKQ/kHRWQHf3vZL2Suvz0FN8HoCCRNXHl1ZWNbtwcOSx/jv2HZ3oPZGdNF0uL0u6zszONzOTdIOkY9ksC0CZRm1MjrvXk5Oe5UlTQ39S0oOSnpL0dP+99ma0LgAlCquPD1tdO6Xb9x0NvSOUSYjlSdXl4u53S7o7o7UAqIjhaYZh9fCBsCFcTEIsD3eKAg0V5xLlOKI2SIex4Zkv7hQFWixO22Fc48ovEhueVUFABxooy3G2w7cORWHDsxoI6EADxR1nG9fO7T0dmp/Tvbdew4ZnhTFtEWigUcfyhw3X2S/udmQmrZxci6y5s+FZbWyKAg0UNtrWJLnWNzAHGfW48bccxa+GIk6KAqiYYMa9sbNBJ06unQ7m0pkN0o2dDSMHZQUvZkb1EdCBhghm5Sura+p2pvTO8zs6cXLtrOeurp0ae9mERPdK3RDQgYaI6myJE7ij0L1SL3S5AA0xaTY93e2M7C+ne6V+yNCBGgo7BRrV2TLd7eiNt94+K1Pvdqb0uVveI0kTdbmg2uhyASoq6uh+WAdLtzOlj1zb00OHl855fM+ubZJoNawzulyAGhoE8aWV1dDOFCm6Vv7Yc8vas2tbZOAmgDcfAR2oiGDmHfy386CNcNTlE/cceJ7su8XYFAUqIizzDhpk3lHSDOFC/RHQgYqI06UyKKOMu3wiyRAu1B8lF6AEk3SpDAzaCONcPsGBoHYiQwcKFjWr/PqrNp2TeVv/595096y5KoPph1EjbTkQ1E5k6EDBknaphLlzx9bQFkYOBLUTAR3IUVhpZdSs8p3bexN1qDDOFsMI6EBOgm2Ig9LKdMiwLCl5mWTSvwTQXAR0ICdRpZVfOW+Dup2picokWV34jGZjUxTISVRp5Wera6fv6DSdu+EZlOWFz2g2MnQgJ6OugZukTDLqwmeydAwjQwdyEnYAKEkHStYXPqO5COhATnZu701UWokStVlKrzmCKLkAORourQw2Nu/Yd3SijU16zREXAR0oQFQLozR+rC295oiLgA4UIO3GJr3miIOADsSUphd81MYmPebICgEdiCFNyUSKbmG8uNtJ9b7AsFRdLmY2bWYPmtlzZnbMzH4/q4UBVTKqZLL/yJJmFw7qyvlHNbtwMPTAT1QLo5ki3xeYVNq2xX+U9A13v0rS70k6ln5JQPWMuvYtzinOqBbGlZCZLqM+DxglccnFzC6S9H5Jfy5J7v6mpDezWRZQnkkun5gyi73ZGbaxGXVJBT3mSCJNhv5bkpYl/bOZHTGzL5rZBRmtCyjFpJdPnPLgVc7rllZWR5ZgBrI6TQpI6QL6eZLeK+mf3H27pF9Kmg8+ycx2m9mimS0uLy+n+Dggf+MunxjcEGSSwkP5GXEGaWV1mhSQJPOIDGPsC81+Q9IT7r6l//s/lDTv7jdFvWZmZsYXFxcTfR6QhXEtglfOPxoaqE3SCwvr/9eeXTg48u7PML3prg7Nz6VYOdrMzA67+8y45yXO0N39fyX92MwG/za8QdKzSd8PyFucMbRx5qYk2bBkkxNFSNuH/leSvmpm75D0I0l/kX5JQHaGM/INZufUvIMbmHHmpkRtkA7KMWxyoiyp2hbd/ai7z7j777r7Tnc/kdXCgLSCGXnUBuZw9hynpj1qI5NNTpSJk6JorLANzjDB7Hnc3JQ4w7I4yo8yENDRWHHq1kmz51FBn0FaKAsXXKCxourWU2a0CKKRyNDRWFEbnARxNBUBHY2VxcUQjLZFnRDQ0Whp6tlpR+YCRaOGDkQYNTIXqCICOhBh1C1DQBUR0IEIccYAAFVCQAcicOoTdcOmKGqlyK6TLLpkgCIlHp+bBONzkUaw60Q6M5e8R0siGizu+FwydNRGWNfJIB2ZtKWQlkQ0ETV01Ma47pJJWgppSUQTEdBRG3G6S+K2FNKSiCYioKM0+48saXbhYKzLlKXwrpOguC2FtCSiiQjoKEWc6+CChi+fkNY3RIdN0lJISyKaiE1RlGJUDXvc5RKDP0/TpUJLIpqIgI5SZFHDTnuRBBdRoGkouaAU1LCB7JGhI3dhpZGoyyeoYQPJkaEjV1Gbn5JOb3ByHRyQDTJ05GKQlS+F1MQHm5+H5ucI4ECGCOjIXNjMlSAO8ADZI6Ajc2EtiUHjNj8ZnAVMjoCOzI3LvsdtfjI4C0iGTVFkblT2HWfzk8FZQDJk6JjYuHJIVEti3C4WBmcByRDQMZE45ZC0x+ovm+6Gdsdw6AgYjYCOicSdwZLkWP1wq+PgJqIBDh0B4xHQMZG8yiHBzN+V7no5oI0I6AgVVSfPqxwSdb1cb7qrQ/Nzqd4baIvUXS5mNmVmR8zskSwWhPKNmlWe1xxxNkKB9LLI0G+TdEzSRRm8FypgVJ18kC2HZe9pDgOxEQqklyqgm9nlkm6S9PeSPpPJilC6cdly2IZn2sNATF8E0ktbcrlX0mclvR31BDPbbWaLZra4vLyc8uNQhCSzytMeBhq+Xo7pi0AyiTN0M7tZ0nF3P2xmH4h6nrvvlbRXkmZmZjzqeaiOJNlyFW4gAtouTYY+K+kWM3tR0gOS5szsvkxWhVIlyZa5gQgon7mnT5r7GfrfuvvNo543MzPji4uLqT8P1RM2MneS4/4AopnZYXefGfc8+tCRibTH/QGkl0mGHhcZOgBMLm6GzvhcAGgISi44jVuCgHojoEMStwQBTUDJBZK4JQhoAgI6JDEcC2gCSi4tN6ibR/U6cTAIqA8CeouFHQYaxnAsoF4I6C0WVjcf4JYgoH4I6C00fHdnGJO4JQioIQJ6y4wrs0jUzYG6osulZUaVWSTq5kCdkaG3zKg2ROrmQL0R0Fsm6u7O3nSXujlQc5RcWubOHVvV7Uyd9RhlFqAZyNBbhrnlQHMR0FuIuzuBZiKgVxjjbAFMghuLSjIuWIf1i5skF90oQNtwp2iFxZk9HtYvPvirl1nlAMLQ5VKCOLPHx42tZVY5gCACegnizB6Pc/yeWeUAhhHQSxAVrIcfD+sXj/s+ANqJgF6COId7dm7vac+uber1g7YF3oPDQACC6HIpyXCXy8XdjsyklZNrke2JtDAC7RW3y4WAnrFJA29Ye2K3M6U9u7YRsAFIih/QKblkaBCcl1ZW5TrTXrj/yFLka+J0vABAHPShJxSWiY8KzlHZdpyOFwCIg4CeQNTBoKiLI0YF56hxtnSwAJgUJZcEojLxKQv2oqwbFZwZZwsgKwT0BKIy7lPu5wRn03oGP7twMLSWPtyeaFqf08KGKIAk6HJJYHbhYOStP4Na+tLK6ulhWgOdDaYLN543sj0RAIJy73IxsyvM7DEzO2Zmz5jZbUnfq25GlUl2bu/p0PycetNdBf+qXHvbdeLkWuwOGACYRJqSy1uS/sbdf1vSdZI+bWa/k82yqi1OmSROlwrtiQCylLjLxd1fk/Ra/9e/MLNjknqSns1obZU2fOvPoIXxjn1HT5dSorpXgmhPBJCVTNoWzWyLpO2Snszi/aoqrPdcUmgL40eu7emhw0uRrYwDtCcCyErqgG5mF0p6SNLt7v7zkD/fLWm3JG3evDntx5Umqvd8Y2dDaAvjY88ta8+ubWfNa/nlm29p7dSZyjrtiQCylKrLxcw6kh6RdMDdPz/u+XXuconqbIlikl5YuOmsxxiwBSCJ3K+gMzOT9CVJx+IE87oaBOFJgrkUXkoZrrsDQNbSdLnMSvqEpDkzO9r/8aGM1lUJw8O2okx3O5z0BFAJabpcHte59y5U2qQlj7Aj/sO6nSl97pb3nH4upRQAZWrNcK6oTU1JE09ClM6cCh28lgAOoGytmeWSZO54VEthb7qrQ/NzBHEAldL4DH3cpmZYFj78muA8FurjAKqq0QE97Hq3oGAWHnyNS6eDerDMAgBV0uiAHmdTM5hth71mEMwPzc/lsUwAyESjA/okm5rjXsPMFQBV1+iAHjUga1S2zZVwAOqq0V0uSa5340o4AHXV6Ax9UE6Z5NBPktcAQBVwBR0AVFzuV9ABAKql0SWXURhlC6BpWhnQk8x1AYCqa2XJJclcFwCoulYGdA4PAWii2pVcsqh9c3gIQBPVKkMfvkHIdab2vf/I0kTvw+EhAE1Uiwx91AjcQe17kiydw0MAmqjyAT3OCNwktW8ubAbQNJUvuYwbgStR+wYAqQYBfVz2Te0bANZVPqCPyr57013t2bWN0gkAqAYBPaoj5d5br+GiZgAYUvlNUTpSACCeygd0iY4UAIij8iUXAEA8BHQAaAgCOgA0BAEdABqCgA4ADVHoJdFmtizppQlecomkn+S0nCrje7cL37tdknzv33T3TeOeVGhAn5SZLca56bpp+N7twvdulzy/NyUXAGgIAjoANETVA/reshdQEr53u/C92yW3713pGjoAIL6qZ+gAgJgqG9DN7EYze97Mfmhm82Wvpwhm9mUzO25mPyh7LUUysyvM7DEzO2Zmz5jZbWWvqQhmttHMvmtm3+9/778re01FMbMpMztiZo+UvZYimdmLZva0mR01s8XM37+KJRczm5L0P5L+SNIrkr4n6WPu/mypC8uZmb1f0uuS/tXdry57PUUxs0slXeruT5nZr0o6LGlnC/73NkkXuPvrZtaR9Lik29z9iZKXljsz+4ykGUkXufvNZa+nKGb2oqQZd8+l/76qGfr7JP3Q3X/k7m9KekDSh0teU+7c/TuSflr2Oorm7q+5+1P9X/9C0jFJjZ+X7Ote7/+20/9RvQwrY2Z2uaSbJH2x7LU0TVUDek/Sj4d+/4pa8B84JDPbImm7pCfLXUkx+qWHo5KOS/qWu7fhe98r6bOS3i57ISVwSd80s8NmtjvrN69qQLeQxxqfubSdmV0o6SFJt7v7z8teTxHc/ZS7XyPpcknvM7NGl9rM7GZJx939cNlrKcmsu79X0p9I+nS/zJqZqgb0VyRdMfT7yyW9WtJaUIB+DfkhSV9194fLXk/R3H1F0n9KurHkpeRtVtIt/VryA5LmzOy+cpdUHHd/tf/zcUn/rvXycmaqGtC/J+ndZnalmb1D0kclfb3kNSEn/c3BL0k65u6fL3s9RTGzTWY23f91V9IHJT1X7qry5e53ufvl7r5F6/9dH3T3j5e8rEKY2QX9TX+Z2QWS/lhSph1tlQzo7v6WpL+UdEDrG2Rfc/dnyl1V/szsfkn/JWmrmb1iZp8qe00FmZX0Ca1na0f7Pz5U9qIKcKmkx8zsv7WexHzL3VvVxtcy75L0uJl9X9J3JT3q7t/I8gMq2bYIAJhcJTN0AMDkCOgA0BAEdABoCAI6ADQEAR0AGoKADgANQUAHgIYgoANAQ/w/7h47kyaTEyQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def f_two_point_line(x1,y1,x2,y2):\n",
    "    '''\n",
    "    两点法确立一个直线，返回所拟合的直线斜率和截距\n",
    "    '''\n",
    "    k = (y2-y1)/(x2-x1)\n",
    "    b = y1 - k*x1\n",
    "    return k,b\n",
    "    \n",
    "\n",
    "def RANSAC_filte_outlier(x,y):\n",
    "    '''\n",
    "    删除离群点\n",
    "    '''\n",
    "    m = 2 #确定一条直线最少每次需要取2个点\n",
    "    e = 0.5 # 只要某个点离所拟合的直线距离小于0.5那么认为不是离群点\n",
    "    n = 20 # 只要某个拟合的直线能够让20个点足够接近它那么认为当前直线是比较理想的直线\n",
    "    K = int(2 * (1/5)**(-m)) # 最大迭代次数，这个可以随便设置的\n",
    "    for _ in range(K):\n",
    "        # 1. 随机选m个点来拟合一条直线这里是两个点所以用两点法确定一条直线\n",
    "        p1_index = np.random.randint(x.shape[0])\n",
    "        p2_index = np.random.randint(x.shape[0])\n",
    "        \n",
    "        # 两点法确定一条直线\n",
    "        p1_x = x[p1_index]\n",
    "        p1_y = y[p1_index]\n",
    "        p2_x = x[p2_index]\n",
    "        p2_y = y[p2_index]\n",
    "        \n",
    "        k,b = f_two_point_line(p1_x,p1_y,p2_x,p2_y)\n",
    "        y_predict = k*x + b\n",
    "        inliner_bool = (np.abs(y_predict-y)<0.5)\n",
    "        if np.sum(inliner_bool)>=20:\n",
    "            return x[inliner_bool],y[inliner_bool]\n",
    "        \n",
    "    return NULL\n",
    "x_filted,y_filted = RANSAC_filte_outlier(x,y)\n",
    "plt.scatter(x_filted,y_filted)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到上图许多噪声点已经被去除。拟合这条曲线可以用最小二乘法求解。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
